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ABSTRACT 

Several properties of the Solar System, including the wide radial spacing and 
orbital eccentricities of giant planets, can be explained if the early Solar System 
evolved through a dynamical instability followed by migration of planets in the 
planetesimal disk. Here we report the results of a statistical study, in which we 
performed nearly 10 4 numerical simulations of planetary instability starting from 
hundreds of different initial conditions. We found that the dynamical evolution 
is typically too violent, if Jupiter and Saturn start in the 3:2 resonance, leading 
to ejection of at least one ice giant from the Solar System. Planet ejection can 
be avoided if the mass of the transplanetary disk of planetesimals was large 
(Mdi S k > 50 Msarth), but we found that a massive disk would lead to excessive 
dynamical damping (e.g., final e 55 < 0.01 compared to present e 55 = 0.044, where 
e 55 is the amplitude of the fifth eccentric mode in the Jupiter's orbit), and to 
smooth migration that violates constraints from the survival of the terrestrial 
planets. Better results were obtained when the Solar System was assumed to 
have five giant planets initially and one ice giant, with the mass comparable to 
that of Uranus and Neptune, was ejected into interstellar space by Jupiter. The 
best results were obtained when the ejected planet was placed into the external 
3:2 or 4:3 resonance with Saturn and M disk ~ 20 M Earth . The range of possible 
outcomes is rather broad in this case, indicating that the present Solar System is 
neither a typical nor expected result for a given initial state, and occurs, in best 
cases, with only a ~5% probability (as defined by the success criteria described 
in the main text). The case with six giant planets shows interesting dynamics 
but does offer significant advantages relative to the five planet case. 
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1. Introduction 

As giant planets radially migrate in the protoplanetary nebula they should commonly 
be drawn into compact systems, in which the pairs of neighbor planets are locked in the 
orbital resonances (Kley 2000, Masset & Snellgrove 2001). The resonant planetary systems 
emerging from the protoplanetary disks can become dynamically unstable after the gas 
disappears, leading to a phase when planets scatter each other. This model can help to 
explain the observed resonant exoplanets (e.g., Gliese 876; Marcy et al. 2001), commonly 
large exoplanet eccentricities (Weidenschilling & Marzari 1996, Rasio & Ford 1996), and 
microlensing data that show evidence for a large number of planets that are free-floating in 
interstellar space (Sumi et al. 2011; but see Veras & Raymond 2012). 

The Solar System, with the widely spaced and nearly circular orbits of the giant planets, 
bears little resemblance to the bulk of known exoplanets. Yet, if our understanding of physics 
of planet-gas-disk interaction is correct, it is likely that the young Solar System followed the 
evolutionary path outlined above. Jupiter and Saturn, for example, were most likely trapped 
in the 3:2 resonance (Masset & Snellgrove 2001, Morbidelli & Crida 2007, Pierens & Nelson 
2008, Pierens & Raymond 2011, Walsh et al. 2011), defined as P Sa t/Pju P = 1-5, where P Jup 
and Psat are the orbital periods of Jupiter and Saturn (this ratio is 2.49 today). 

To stretch to its current configuration, the outer Solar System most likely underwent 
a dynamical instability during which Uranus and Neptune were scattered off of Jupiter and 
Saturn and acquired eccentric orbits (Thommes et al. 1999, Tsiganis et al. 2005, Morbidelli 
et al. 2007, Levison et al. 2011). The orbits were subsequently stabilized (and circularized) 
by damping the excess orbital energy into a massive disk of planetesimals located beyond 
the orbit of the outermost planet (hereafter transplanetary disk), whose remains survived 
to this time in the Kuiper belt. Finally, as evidenced by dynamical structures observed in 
the Kuiper belt, planets radially migrated to their current orbits by scattering planetesimals 
(Malhotra et al. 1995, Gomes et al. 2004, Levison et al. 2008). 

Several versions of the Solar System instability were proposed. Thommes et al. (1999), 
motivated by the apparent inability of the existing formation models to accrete Uranus 
and Neptune at their present locations, proposed that Uranus and Neptune formed in the 
Jupiter-Saturn zone, and were scattered outwards when Jupiter, and perhaps Saturn, ac- 
creted nebular gas. This work represented an important paradigm shift in studies of the 
Solar System formation. Inspired by this work, Levison et al. (2001) suggested that the 
instability and subsequent dispersal of the planetesimal disk, if appropriately delayed (or 
due to the late formation of Uranus and Neptune; Wetherill 1975), could explain the Late 
Heavy Bombardment (LHB) of the Moon. The LHB was a spike in the bombardment rate 
some 4 Gyr ago suggested by the clustering of ages of several lunar basins. The nature of the 
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LHB, however, is still being debated (see Hartmann et al. 2000 and Chapman et al. 2007 
for reviews). 

In the Nice version of the instability (Tsiganis et al. 2005), Uranus and Neptune were 
assumed to have formed just outside the orbit of Saturn, while Jupiter and Saturn were 
assumed to have initial orbits (where "initial" means at the time when the protoplanetary 
nebula dispersed), such that Psat/-Pjup < 2. The instability was triggered when Jupiter 
and Saturn, radially drifting by scattering planetesimals, approached and crossed the 2:1 
resonance (-Psat/-Pj up = 2). The subsequent evolution was similar to that found in the 
Thommes et al. model, but led to a better final orbital configuration of the planets (assuming 
that the disk contained ~ 20-50 Msarth and was truncated at 30-35 AU). 

The Nice model is compelling because it explains the separations, eccentricities, and 
inclinations of the outer planets (Tsiganis et al. 2005), and many properties of the popula- 
tions of small bodies (Morbidelli et al. 2005, Nesvorny et al. 2007, Levison et al. 2008, 2009, 
Nesvorny & Vokrouhlicky 2009). It is currently the only migration model that is consistent 
with the current dynamical structure of the terrestrial planets (Brasser et al. 2009) and the 
main asteroid belt (Morbidelli et al. 2010). Moreover, the Nice model can also reproduce 
the magnitude and duration of the LHB (Gomes et al. 2005, Bottke et al. 2012). 

Two potentially problematic issues of the original Nice model (hereafter ONM) were 
addressed by Morbidelli et al. (2007) and Levison et al. (2011). Morbidelli et al. pointed 
out that the initial planetary orbits in the ONM were chosen without the proper regard to the 
previous stage of evolution during which the giant planets interacted with the protoplanetary 
nebula. As shown by hydrodynamical studies (Masset & Snellgrove 2001, Morbidelli & Crida 
2007, Pierens & Nelson 2008, Pierens & Raymond 2011), this interaction most likely led to 
the convergent migration of planets, and their trapping in orbital resonances. Morbidelli et 
al. (2007) studied the dynamical instability starting from the resonant planetary systems 
and showed that the orbit evolution of planets was similar to that of the ONM. 

To produce the LHB, the instability needs to occur late relative to the dispersal of the 
protoplanetary nebula. A protoplanetary nebula typically disperses in ~3-10 Myr after the 
birth of the star (Haisch et al. 2001), which for the Sun was probably contemporary to the 
formation of the first Solar System solids, 4.568 Gyr ago (Bouvier et al. 2007, Burkhardt 
et al. 2008). The onset of the LHB, traditionally considered to be 3.9-4.0 Gyr ago (Ryder 
2002), has been recently revised to 4.1-4.2 Gyr ago (Bottke et al. 2012). This means that 
the instability had to occur with a delay of 350-650 Myr, with the exact value depending on 
the actual time of the LHB event Such a late instability can be triggered in the ONM if 



Here we opted for including the full range of previous and new estimates of the delay, because it is not 
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the inner edge of the planetesimal disk was close, but not too close, to the outer ice giant. If 
the edge had been too close, the instability would have happened too early to be related to 
the LHB. If the edge had been too far, the instability would have happened too late or would 
not have happened at all, because the planetesimals would had stayed radially confined and 
not evolved onto planet-crossing orbits (where they could influence planets by short-range 
interactions). 

To avoid the need for fine tuning of the inner disk's edge in the ONM, Levison et al. 
(2011) proposed that the late instability was caused by the long-range interactions between 
planets and the radially-confined distant planetesimal disk. The long-term interactions arise 
when gravitational scattering between disk's planetesimals is taken into account. As a result 
of these interactions, the planets and planetesimal disk slowly exchanged the energy and 
angular momentum until, after hundreds of Myr of small changes, the resonant locks between 
planets were broken and the instability reigned supreme. Note that the instability trigger 
in Levison et al. (broken resonant locks) is fundamentally different from that of the ONM 
(major resonance crossing during migration). 

Here we report the results of a new statistical study of the Solar System instability. 
Sections 2-4 explain our method and constraints. The results are presented in Sect. 5. We 
discuss the plausible initial configurations of planets after the gas nebula dispersal, mass of 
the planetesimal disk, and effect of different trigger mechanisms on the results. The first 
steps in this direction were taken by Batygin & Brown (2010), who found that some initial 
resonant configurations may work while others probably do not. 

Following Nesvorny (2011) and Batygin et al. (2012) we also consider cases in which 
the young Solar System had extra ice giants initially and lost them during the scattering 
phase (these works and their relation to the present paper will be discussed in more detail 
below; see, e.g., Sect. 6). We extend the previous studies of the Solar System instability by: 
(1) exploring a wide range of initial parameters (new resonant chains, six planets, etc.), (2) 
improving the statistics with up to 100 simulations per case, (3) considering different trigger 
mechanisms (including the one recently proposed by Levison et al. 2011), and (4) applying 
strict criteria of success/failure to all studied cases (see Sect. 3 and 4). Our conclusions are 
given in Sect. 7. 



100% guaranteed that the new estimate is correct. Imposing the 350-450 My delay, instead of 350-650 My, 
would not make much of a difference, because most planetary systems that are stable over 450 My are also 
stable over 650 My. 
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2. Method 

We conducted computer simulations of the early evolution of the Solar System. In 
the first step (Phase 1), we performed hydrodynamic and iV-body simulations to identify 
the resonant configurations that may have occurred among the young Solar System's giant 
planets. Our hydrodynamic simulations used Fargo (Masset 2000) and followed the method 
described in Morbidelli et al. (2007). As the Fargo simulations are CPU expensive, we used 
these results as a guide, and generated many additional resonant systems with the iV-body 
integrator known as SyMBA (Duncan et al. 1998). 

Planets with masses corresponding to those of Jupiter, Saturn and ice giants, ordered in 
increasing orbital distance from the Sun, were placed on initial orbits with the period ratios 
slightly larger that those of the selected resonances. We tested several cases for the initial 
radial order of ice giants including the case with Uranus on the inside of Neptune's orbit, 
Uranus on the outside of the Neptune's orbit, and case where both ice giants were given 
masses intermediate between that of Uranus and Neptune. The results (discussed in Sect. 
5) obtained in these cases were statistically equivalent showing that the initial ordering of 
Uranus and Neptune was not important. 

The planets were then migrated into resonances with SyMBA modified to include forces 
that mimic the effects of gas. We considered cases with four, five and six initial planets, 
where the additional planets were placed onto resonant orbits between Saturn and the inner 
ice giant, or beyond the orbit of the outer ice giant. Additional planets were given the mass 
between 1/3 and 3 times the mass of Uranus. 

Different starting positions of planets, rates of the semimajor axis and eccentricity evo- 
lution (as implied by different gas disk densities), and timescales for the gas disk's dispersal 
produced different results. For Jupiter and Saturn, we confined the scope of this study to 
the 3:2 and 2:1 resonances, because the former one is strongly preferred from previous hy- 
drodynamic studies (Masset & Snellgrove 2001, Morbidelli et al. 2007, Pierens & Nelson 
2008, Walsh et al. 2011, Pierens & Raymond 2011). The 2:1 resonance was included for 
comparison. According to our tests with Fargo, trapping of Jupiter and Saturn in the 2:1 
resonance may require special conditions (e.g., a low gas density implying slow convergent 
migration; Thommes et al. 2008). The 2:1 resonance does not generically lead to the out- 
ward migration of Jupiter and Saturn that is required for the Grand Tack model (Walsh et 
al. 2011). 

The planetary eccentricities, inclinations and resonant amplitudes obtained at the end 
of our Phase-1 simulations were e < 0.1, % < 0.2° and < 60°. The inner ice giant typically had 
the most eccentric orbit (0.05 < e < 0.1), while the other planets' orbits were more circular 
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(e.g., Jupiter had e < 0.01 in most cases). We also considered cases with several different 
damping strengths for the same resonant chain. These cases helped us to understand the 
effect of the initial excitation on the evolution of orbits during the instability. Changing the 
damping strength may not be strictly justified, because the semimajor axis migration and 
eccentricity damping are expected to be coupled (Goldreich & Sari 2003). The exact nature 
of this coupling, however, depends on several disk parameters that are poorly constrained 
(e.g., Lee & Peale 2002). 

The instability of a planetary system can occur after the gas disk's dispersal when the 
stabilizing effects of gas are removed. Such an instability can be triggered spontaneously 
(e.g., Weidenschilling & Marzari 1996), by divergent migration of planets produced by their 
interaction with planetesimals leaking into the planet-crossing orbits from a transplanetary 
disk (Tsiganis et al. 2005), or by long-range perturbations from a distant planetesimal disk 
(Levison et al. 2011). 

It is often assumed that the instability in the Solar System occurred at the time of 
the LHB. As discussed in the introduction this requires that the giant planets remained on 
their initial resonant orbits for ~350-650 Myr. To allow for this possibility, we examined the 
resonant configurations identified above and selected those that were stable over hundreds of 
Myr, if considered in isolation (i.e., no disk, planets only, no external perturbations). Only 
the stable systems were used for the follow-up simulations, in which we tracked the evolution 
of planetary orbits through and past the instability (Phase 2). 

We included the effects of the transplanetary planetesimal disk in the Phase-2 simula- 
tions. The disk was represented by 1,000 equal-mass bodied that were placed into orbits with 
low orbital eccentricities and inclinations (0.01), and radial distances between r in < r < r out . 
We also performed a limited number of simulations with 100, 300, 3,000, and 10,000 disk 
bodies to test the effect of resolution^ and with excited disks, where the eccentricities and 
inclinations of disk bodies were set to be ~0.1 (Sect. 5.2.1). The surface density was fixed 
at S(r) = 1/r in all simulations, because our tests showed that considering different density 
profiles has only a minimal effect on the results (see also Batygin & Brown 2010). The outer 



2 Planets fully interact with each other and with the disk bodies. The gravitational interaction between 
disk bodies was neglected to cut down the CPU cost. This is justified because the behavior of the sys- 
tem during the instability is mainly driven by the interaction between planets, and between planets and 
planetesimals. 

3 Note that even our highest resolution runs did not have enough bodies to properly model the planetesimal 
disk, which is thought to have contained ^1,000 Pluto-size and myriads of smaller planetesimals (Morbidelli 
et al. 2009a). This should not be a major problem, however, because our results with increased resolution 
showed only a minor dependence on the number of bodies used. 
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edge of the disk was placed at r out = 30 AU, so that the planetesimal-driven migration is 
expected to park Neptune near its present semimajor axis (Gomes et al. 2004)0 

We considered cases with r in = a n + A, where a n is the semimajor axis of the outer ice 
giant, n is the number of planets, and A = 0.5-5 AU. The instability was usually triggered 
early for A < 1 AU. To trigger the instability for A > 1 AU, we broke the resonant locks by 
altering the mean anomaly of one of the ice giants. This was done at the start of Phase-2 
simulations. This method was inspired by the results of Levison et al. (2011). Different 
cases are denoted by B in the following text, where B(j) stands for breaking the resonant 
lock of ice giant planet j (e.g., B(l) corresponds to the innermost ice giant; £>(0) is the ONM 
trigger) . 

In either case discussed above, the scattering phase between planets started shortly 
after the beginning of our simulations, which guaranteed low CPU cost. We considered 
different masses of the planetesimal disk, Mdi S k, with Mdi S k between 10 and 100 Msarth- 
Thirty simulations were performed in each case, where different evolution histories were 
generated by randomly seeding the initial orbit distribution of planetesimals. The number 
of simulations was increased to 100 in the interesting cases. In total, we completed nearly 
10 4 scattering simulations from over 200 different initial states. The new trials were selected 
with the knowledge of the results of the previous simulations and were optimized to sample 
the interesting parts of parameter space. Each system was followed for 100 Myr with the 
standard SyMBA integrator (Duncan et al. 1998), at which point the planetesimal disk was 
largely depleted and planetary migration ceased. We used a h = 0.5 yr time step and verified 
that the results were statistically the same with h = 0.25 yr. 



3. Constraints 

We defined four criteria to measure the overall success of simulations. First of all, the 
final planetary system must have four giant planets (criterion A) with orbits that resemble 
the present ones (criterion B). Note that A means that one and two planets must be ejected in 
the five- and six-planet planet cases, while all four planets need to survive in the four-planet 
case. As for B, we claim success if the final mean semimajor axis of each planet is within 20% 
to its present value, and if the final mean eccentricities and mean inclinations are no larger 
than 0.11 and 2°, respectively. These thresholds were obtained by doubling the current mean 



4 The real planetesimal disk probably continued all the way to ~47 AU, as implied by the existence of the 
cold classical Kuiper belt (Batygin et al. 2011), but this extension likely contained too little mass to drive 
Neptune's migration. 
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eccentricity of Saturn (es a t = 0.054) and mean inclination of Uranus (iura = 1-02°; Table 1). 

For the successful runs, as defined above, we also checked on the history of encounters 
between giant planets, evolution of the secular g 5 , g 6 and s 6 modes, and secular structure of 
the final planetary systems. To explain the observed populations of the irregular moons that 
are roughly similar at each planet (when corrected for observational incompleteness; Jewitt 
& Haghighipour 2007), all planets -including Jupiter- must participate in encounters with 
other planets (Nesvorny et al. 2007). Encounters of Jupiter and/or Saturn with ice giants 
may also be needed to excite the g§ mode in the Jupiter's orbit to its current amplitude 
(ess = 0.044; Morbidelli et al. 2009b). Tables 2-4 list the secular frequencies and amplitudes 
of the outer Solar System planets. 

It turns out that it is generally easy to have encounters of one of the ice giants to 
Jupiter if the planets start in a compact resonant configuration (such encounters occur in 
most simulations for most cases tested here). The amplitudes of g% and sq modes also do not 
pose a problem. It is much harder to excite e^, however. We therefore opt for a restrictive 
criterion in which we require that e 55 > 0.022 in the final systems, i.e., at least half of 
its current value (criterion C). The e 55 amplitude was determined by following the final 
planetary systems for additional 10 Myr (without planetesimals), and Fourier- analyzing the 
results (Sidlichovsky & Nesvorny 1996). 

The evolution of secular modes, mainly g 5 , g 6 and s 6 , is constrained from their effects 
on the terrestrial planets and asteroid belt. As giant planets scatter and migrate, these 
frequencies change. This may become a problem, if g 5 slowly swipes over the g\ or g 2 modes, 
because the strong g\ = g$ and g 2 = g§ resonances can produce excessive excitation and 
instabilities in the terrestrial planet system (Brasser et al. 2009, Agnor & Lin 2011). The 
behavior of the g§ and s$ modes, on the other hand, is important for the asteroid belt 
(Morbidelli et al. 2010, Minton & Malhotra 2011). 

As (?5, <?6 and sq are mainly a function of the orbital separation between Jupiter and 
Saturn, the constraints from the terrestrial planets and asteroid belt can be conveniently 
defined in terms of -Psat/-Pjup- This ratio needs to evolve from <2.1 to >2.3 in < 1 Myr 
(criterion D; also final -Psat/-Fjup < 2.8; i.e., close to present -Psat/-Pj up = 2.49 but not too 
restrictive; see Sect. 4), which can be achieved, for example, if planetary encounters with 
an ice giant scatter Jupiter inward and Saturn outward. The 1-Myr limit is conservative in 
that a slower evolution of secular frequencies, which fails to satisfy the criterion D, would 
clearly violate the constraints. Note also that in most simulations that satisfy criterion D, 
Psat/-Pjup evolves from <2.1 to >2.3 in much less than 1 Myr. Most simulations that satisfy 
D should therefore satisfy the constraints from the terrestrial planets and asteroids. 
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In summary, our constraints A and B express the basic requirements on a successful 
model. Constraint C is more restrictive in that it places a more precise condition on the 
dynamical structure of the final systems. Constraint D is the least rigorous. Given that 
it is related to the orbits of terrestrial planets, this constraint would not need to apply if 
the giant planet instability occurred before the terrestrial planets formed (i.e., early and well 
before the LHB). Note, however, that evolutions violating the constraint D would be difficult 
to reconcile with the present dynamical structure of the asteroid belt (Walsh & Morbidelli 
2011), regardless the timing of the giant planet instability. We therefore do not see a way 
around criterion D. 

We report our results in Sect. 5. In doing so, we consider criteria C and D independently, 
but also discuss the most interesting cases where both C and D were satisfied simultaneously. 
Note that the success rates for C, D and C&D were computed over the subset of simulations 
that satisfied A and B (i.e., it would not make sense, for example, to claim success for C 
if the system ended with incorrect number of planets). Also, B can be satisfied only if 
A is satisfied. The reported fractions were normalized by the total number of simulations 
performed in each case (equal to 30 or 100). 

We ignored other constraints from the populations of small bodies in this work (e.g., 
cold classical Kuiper belt survival; Batygin et al. 2012). These constraints will be addressed 
in the upcoming work. 

4. Measure of Success and Failure 

The instability is a chaotic process and as such it is not expected to produce deterministic 
results. Two planetary systems starting from the same initial conditions (that is, initial after 
the dispersal of the protoplanetary nebula) are expected to follow different evolutionary 
paths during the instability, and end up producing different final systems. It would thus 
be worthless to base our conclusions on one or a small number of simulations that may 
not be representative of the full range of possible outcomes. Here we performed 30 to 100 
simulations for each initial setup. The number of simulations was chosen as a compromise 
between the need for reasonable statistics and our CPU limitations. 

A fundamental difficulty with studying the Solar System instability is that the Solar 
System, as we know it now, is just one realization of all possible evolution paths. We do not 
know if this realization is typical of those arising from the initial state, or the nature played 
a trick on us, and the Solar System has taken an unusual path. The latter possibility would 
offer a grim perspective on the possibility to determine the Solar System's state before the 
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instability. Here we assume that the Solar System followed a path that was not "exceptional" . 
Therefore we assume some cutoff probability for success in matching the criteria described 
in Sect. 3, and give preference to those initial conditions and evolution paths that exceed 
the cutoff. Now, the assumed cutoff must depend on how selective our criteria are because it 
is obvious that the measure of the set of initial conditions that lead to the exact architecture 
of the Solar System is practically zero. We therefore cannot set the criteria that would be 
too restrictive because we would not be able to obtain satisfactory statistics with our 30-100 
simulations per case. The criteria described in Sect. 3 were defined with this in mind. 

For example, there is nothing wrong with the initial setup that would lead to 3 or 
4 giant planets in the end, as long as the two results occur with a roughly comparable 
likelihood. Statistical fluctuations are also expected for criteria B and C, because some 
systems may become more excited and/or radially spread (thus failing on B), while others 
end up dynamically cold (leaving Jupiter with low and failing C). Criterion D can be 
even more difficult to satisfy because planetary encounters produce a rather large variety 
of jumping Jupiter histories, in some of which Jupiter does not jump enough (leading to 
-Psat/-Pj up < 2.3) or jumps too much (leading to -Ps a t/-Pju P 3> 2.5). Both these cases would 
fail D. 

If we would optimistically assume that the success rate of 50% for each of our four 
criteria would be satisfactory, the combined probability to simultaneously satisfy all of them 
would be 1/16 ~ 6%. Based on this we set the target cutoff probability to be of the order 
of a few percent. The number of simulations used in this work (see Sect. 2) was chosen so 
that we are able to measure such probabilities. 

Note that in many cases studied here the differences between the model results and 
present Solar System are systematic. For example, in the runs starting with four planets, 
e 55 systematically ends up too low, because the evolution needs to be relatively tranquil to 
satisfy A and B, and fails on C because is not excited enough (Sect. 5.1.1). It is clear 
that the probability to match our criteria is 1% in these cases (although we do not have 
enough statistics to say how low the success rate actually is). 

5. Results 

5.1. Case with 4 Planets 

We start by discussing the results of simulations with four initial giant planets. All 
four planets have to survive in this case. We performed 2670 integrations for many different 
resonant chains, including cases of the 2:1, 3:2, 4:3 and 5:4 resonances between pairs of 
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planets. Table 5 summarizes the results of selected simulations. These cases were the most 
interesting ones, i.e. those where the criteria defined in Sect. 3 were satisfied with the highest 
percentages (Tables 6 and 7 show the selected simulations for the 5- and 6-planet cases). 
The second and higher order resonances were not considered. The results can be divided 
into two broad categories that share common traits: (1) Jupiter and Saturn starting in 3:2, 
and the inner ice giant in the 3:2 or higher degree resonance with Saturn (e.g., 4:3); and (2) 
Jupiter and Saturn in 3:2, and the inner ice giant in the 2:1 resonance with Saturn (or any 
other comparably distant orbit), or Jupiter and Saturn in the 2:1 resonance. We discuss the 
two categories in sections 5.1.1 and 5.1.2, respectively. 



5.1.1. 3:2 and Tight 

In category (1), the system can easily be destabilized, and when it becomes unstable, the 
instability tends to be violent: one ice giant gets ejected from the Solar System or both ice 
giants survive, but their orbits end up being very excited (e > 0.1) and/or scattered to large 
heliocentric distances (a > 50 AU). This happens particularly for disk masses M d i S k ^$ 35 
MEarth- Thus the simulations with low disk masses are clearly unsuccessful in matching the 
present Solar System. 

For example, in the case with (3:2,3:2,4:3) very closely matching one of the systems 
discussed in Morbidelli et al. (2007), 87% of simulations show ejection of one or more ice 
giants (for M disk = 35 M^ai-th and £>(1)). None of the remaining 13% satisfied our criterion B, 
typically because the final orbits were overly excited (e.g., Uranus ended up with i ~ 10°). 
The simulations with M disk < 35 M Earth show even larger percentage of ejected planets and 
more excitation than the ones with M disk = 35 M Earth . Thus, disks with M disk < 35 M Earth 
do not apply. This result is robust in that it does not depend on the initial resonant chain 
(within the definition of category (1)). 

A single exception from the above rule was one simulation started with (3:2,4:3,3:2), 
which was, perhaps incidentally, a resonant chain favored by Batygin & Brown (2010). In 
this particular run, Neptune was scattered to a ~ 26 AU, where it was quickly circularized 
by dynamical friction from the planetesimal disk (Fig. [1]). The subsequent migration of 
Neptune and Uranus was minimal and the two orbits ended too close to the Sun (a = 16 
and 26 AU, respectively). 

On the positive side, Jupiter's orbit was excited by the encounters with Uranus, such 
that ej up = 0.042 and 2j up = 0.37° in the end (bar indicates mean values; initial values 
were 0.002 and 0.013°, respectively), both in a good agreement with the present values 
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(ejup = 0.046 and zj up = 0.37°; Table 1). Even more encouragingly, final = 0.036, also 
in a close agreement with present = 0.044 (Table 3). The main problem with this run, 
however, was that the jump in Jupiter and Saturn orbits produced by scattering encounters 
with Uranus and Neptune was minimal (Psat/-fjup moved to ~ 1.8 at t < 10 5 yr) so that the 
final -Psat/-Pju P ratio was too low (~2.1). 

It is difficult to evaluate whether this is a fundamental problem, or whether we would 
find better cases if we increased the number of simulations. We increased the number of 
simulations to 100 but did not obtain a case that would be even remotely as good as the 
one shown in Fig. [TJ We were therefore probably just lucky when finding this case with 
fewer runs. We conclude that it might be possible to match the A, B and C constraints with 
Miisk ~ 35 M Earth , but the likelihood of this occurring is <1%. If, in addition, criterion D is 
considered, the probability drops to <Cl%. 

A better success in matching criteria A and B was obtained with more massive disks 
(Mdisk > 50 MEarth)- The success for A and B increased with M^sk, because the more massive 
disks were capable of stabilizing the system of four planets more efficiently. The best success 
rate for A and B was obtained with M disk = 75 M Earth and r out = 26 AU, where 67% of 
jobs matched A and 40% satisfied B. We set the outer disk edge at 26 AU, because in our 
initial runs with r out = 30 AU Neptune migrated to ~35 AU. A massive planetesimal disk 
is apparently capable of migrating Neptune beyond the initial outer disk's edge. 

The overall distribution of final orbits was reasonable in this case (Fig. [2]). Uranus 
ended up slightly inside its present orbit (aura = 17.5 AU compared to real au ra = 19.2 AU), 
which may or may not be a problem. On one hand, Uranus ended up too close to the Sun 
in most of the simulations that we performed with 4 planets and the 3:2 resonance between 
Jupiter and Saturn. So, this is a systematic effect. On the other hand, one might be able to 
resolve this problem by modifying the initial setup (e.g., starting Uranus on a larger orbit; 
but see Sect. 5.1.2). Alternatively, the present orbit of Uranus is an outlier in the semimajor 
axis distribution about 2a from the mean. 

A more fundamental (and also systematic) problem of these simulations was that both 
Jupiter and Saturn ended up on orbits that were too circular. The mean eccentricities that 
we obtained in the simulations with Mdi S k = 75 MEarth were ej up = 0.01 and es a t = 0.02, 
while real ej up = 0.046 and es a t = 0.054. This was a consequence of the damping effect 
of massive planetesimal disk on the gas giant orbits (see Fig. [3] for an illustration). This 
problem was general for all simulations that we performed with massive disks. 

To understand this effect better, we analyzed the simulations with massive planetesimal 
disks in more detail. We found that it is generally not a problem to excite ej up and es a t to 
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~0.05 by scattering encounters with one of the ice giants. What is problematic, however, is 
to maintain the excited state when the massive planetesimal disk is present. In fact, in all 
our simulations with massive disks, ej up and es a t were quickly damped after the excitation 
event. 

We believe that the damping effect arises from secular friction (Levison et al. 2011) and 
occurs when e and/or the longitude of perihelion w of a planet suddenly changes (e.g., by 
being scattered off of another planet or due to orbital resonance crossing). This changes the 
phase and amplitude of the secular forcing between the planet and the planetesimal disk. 
When the system relaxes the mean eccentricity of the planetesimals increases, while that of 
the planet decreases. The secular friction damps planet's inclination as welljj 

The excessive damping of ej up and e$ a t could potentially be avoided if the planetesimal 
disk were dynamically excited or depleted prior to the event that excites 655. This might 
happen, for example, if Neptune were scattered deep into the planetesimal disk (to 28-30 
AU) and disrupted it before the ess excitation event. Unfortunately, we did not see this 
happening in our simulations, except for a few cases that were flawed for other reasons (e.g., 
Fig. [T]). Another possibility, discussed in more detail in Sect. 5.2.1, is that the planetesimal 
disk was stirred by large objects that formed in the disk and/or by resonances with ice 
giants. Whatever the stirring mechanism was, however, we found that a disk that did not 
exert a strong damping effect on Jupiter and Saturn was not capable of preventing ejection 
of Uranus or Neptune. These simulations thus failed on A or C. 

There were two additional problems with the simulations such as the one illustrated in 
Fig. [3j First, the ones that matched criteria A and B were typically those that avoided any 
strong interactions between planets. This resulted in a nearly smooth migration histories of 
planets that violated D. In the specific (and representative) case showed in Fig. EJ the system 
spent more than 3 Myr with 2.1 < Psat/-Pjup < 2.3, which is expected to be incompatible 
with the survival of the terrestrial planets (as discussed in Sect. 3). 

Second, the smooth crossing of the 2:1 resonance between Jupiter and Saturn excited 
but not 655. The simulations therefore failed on C. The excitation by resonant crossing 
mainly occurs as the eccentricity of Saturn evolves over the resonance with Jupiter (Jupiter's 
eccentricity changes less because Jupiter has a larger mass). Thus, the resonant crossing di- 
rectly affects the amplitudes related to the g§ mode, and not those of the g$ mode (Morbidelli 
et al. 2009b). The needed excitation of e 55 could more easily be achieved by scattering en- 
counters of an ice giant with Jupiter, but such encounters generally lead to violent orbital 



5 Secular friction is different from an effect known as the dynamical friction (Binney & Tremaine 1987), 
which arises from encounters between bodies. 



-14- 



histories that fail A or B if only four giant planets are considered. 

Finally, we discuss the most promising four-planet case obtained with (3:2,3:2,4:3) and 
Mdisk = 50 M Earth . In this case, 27% of runs satisfied A, 11% of runs satisfied B, and we also 
found one case (out of a hundred) that satisfied criterion C (Fig. @J. This specific simulation 
ended with g 5 = 5.98 arcsec yr -1 and g% = 32.4 arcsec yr -1 , reflecting the fact that Jupiter 
and Saturn finished a bit closer to each other than in reality, and = 0.026, = 0.010, 
e 66 = 0.031, and e 65 = 0.022. This is a pretty good match to the real values (Table 3), 
although ess is only a bit larger than a half of the present value. Unfortunately, although 
this simulations also showed additional encouraging signs (e.g., a ~0.5 jump of Ps&t/Pjup), 
it violated D, because the jump of -Psat/-fj up wa s not large enough. 

In summary, we were not able to find any initial condition for the four-planet case, 
where the likelihood of simultaneously matching all our four criteria would exceed 1%. 

5.1.2. 3:2 and Loose, or 2:1 

Here we first consider cases where Jupiter and Saturn are in the 3:2 resonance, the 
inner ice giant in the 2:1 resonance with Saturn, and various resonances between the first 
and second ice giants. These cases showed a different behavior than the ones discussed in 
the previous section. They were more stable in that breaking the resonant locks of the inner 
ice giant did not always generate an instability. If it did, the instability was usually mild 
such that close encounters between planets did not occur, and planets ended up migrating 
smoothly. This had two conflicting consequences. The smooth migration gave a large success 
rate for A and B, because the system did not suffer destabilizing perturbations and evolved 
quietly. In most other aspects, however, these simulations produced incorrect results. 

First of all, Jupiter and Saturn did not reach the current -Ps a t/-Pj up = 2.5 unless the 
disk was heavy. This was because, in absence of encounters with ice giants, the gas giants 
did not jump, and all the migration work had to be accomplished by planetesimals. With 
heavy disks, which gave a more reasonable final -Psat/-Pj up ratio, Psat/-Pj up slowly migrated 
over 2.1-2.3 and violated our constraint D. 

The eccentricities of Jupiter and Saturn were excited when these planets crossed their 
mutual 2:1 resonance. The resonant crossing led to ej up ~ 0.04 and es a t — 0.08, which was 
lovely, except that this eccentricity excitation was contained in the secular modes related to 
g§. The simulations therefore violated C as was never large enough. We have not found a 
single case where the constraint C or D were satisfied. We conclude that the resonant chains 
with the inner ice giant in the 2:1 resonance with Saturn can be ruled out. 
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The case with the 2:1 resonance between Jupiter and Saturn showed behavior that was 
reminiscent to that of the ONM (Tsiganis et al. 2005). The success rate for A and B was 
large, but Jupiter and Saturn migrated smoothly and criterion D was not satisfied. Also, 
because of the lack of an adequate excitation agent, Jupiter ended up with orbital eccentricity 
that was too low (Fig. [5]), and C was not satisfied as well. In addition, Uranus typically 
reached a\j TSu < 18 AU, well inside its present orbit. In the ONM (Tsiganis et al. 2005), the 
inner ice giant was started far beyond the 2:1 resonance with Saturn, which helped Uranus 
to reach 19.2 AU. 



5.2. Case with 5 Planets 

Following Nesvorny (2011) and Batygin et al. (2012) we considered the case with five 
giant planets as an interesting solution to some of the problems discussed in the previous 
section. The main advantage of the five planet case is that the system can afford to lose 
a planet. This resolves, to a degree, the conflicting requirements for a significant excita- 
tion of ess, which requires strong scattering encounters between Jupiter and ice giants, and 
constraints A and B. 

We performed 4050 integrations with five planets for 19 different resonant chains. We 
tested different masses of the fifth planet and different initial orbits. We found that the best 
results were obtained when the fifth planet had mass comparable to that of Uranus/Neptune, 
and was placed between the orbits of Saturn and Uranus. We discuss on this case below. The 
case in which the fifth planet was given a large mass typically resulted in a violent instability 
and ejection of Uranus and/or Neptune from the system. If the fifth planet was given a lower 
mass and orbit beyond that of Neptune, the inner ice giant got ejected during the instability 
thus leaving a final system with incorrect masses. This later case bears similarities to the 
four-planet case discussed in Sect. 5.1. 

5.2.1. (3:2,3:2,4:3,5:4) and Alike 

Morbidelli et al. (2007) found one system with five planets that remained dynamically 
stable over 1 Gyr (in absence of perturbations from the planetesimal disk)@ The five planets 



6 In the published article, Morbidelli et al. (2007) claimed that the five-planet configuration was unstable. 
Unfortunately, to test the stability, Morbidelli et al. (2007) used the original version of SyMBA (Duncan et 
al. 1998), which had a problem with integrating very compact planetary configurations for long periods of 
time. After fixing this problem, Levison et al. (2011) found that the five-planet configuration identified by 
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were in the (3:2,3:2,4:3,5:4) resonant chain. The extra ice giant had mass similar to that of 
Uranus/Neptune. The orbital eccentricities were 0.004, 0.02, 0.08, 0.035, and 0.01 (from the 
inner to outer planet). Levison et al. (2011) showed that in such a case, where the inner ice 
giant had the largest eccentricity, it can be expected that the (late) instability was triggered 
by the inner ice giant (i.e., 23(1) in notation introduced in Sect. 2). We discuss this case first 
and get to different trigger mechanisms and resonant chains later. 

We found that the low-mass disks with M disk < 20 M Earth did not work, because the 
system frequently lost two or more planets, ending with three or less. Specifically, only ~10% 
of the simulations with M disk = 20 MEarth ended up with four planets, but the orbits were 
all wrong. This is similar to what happened in the four-planet case described in Sect. 5.1.1 
for M disk < 35 M Earth . The instability is simply too violent in this case to be contained by a 
low-mass planetesimal disk. Very massive planetesimal disks with M disk > 50 M Earth did not 
work as well mainly because too much mass was processed through the Jupiter/Saturn region, 
leading to excessive migration of Jupiter and Saturn, and incorrect final Psat/-Pj up > 3. 

Intermediate disk masses, M disk = 35 M Earth and 50 M Earth , were more promising. With 
A = 1 AU (recall that A denotes the initial radial separation between the outer ice giant 
and inner edge of the planetesimal disk; see Sect. 2), these cases showed the 13% and 37% 
success rates for A, and 3% and 23% for B, respectively. These fractions are pretty much 
independent of the instability trigger. The larger success rate for M disk = 50 M Eart h than 
for M disk = 35 M Eart h reflects the stabilizing effect of the planetesimal disk, which increases 
with its mass. The success rate tends to drop with increasing A, but this trend is not always 
clear. For example, the success rate for A drops from 37% to 23% when A is increased to 
3.5 AU for Mdisk = 50 M Eartri (similarly, B drops from 23% to 10%), while the percentages 
for A = 1 AU and 3.5 AU are similar for M disk = 35 M Eart h. 

Figure E] shows the final systems for M disk = 50 M Earth . This result is encouraging. 
The distributions show a larger range of values than what we got in the four-planet case, 
because the interaction of five planets is more complex and leads to a larger variety of results. 
Neptune's model orbit tends to be slightly more excited than the real one, but this difference 
is small and probably not fundamental in that it could be resolved by tuning the parameters. 
On the positive side, some of the simulations that satisfied A and B are now also showing 
signs of success for the criterion D (Fig. [7]). This has not happened in the four-planet case 
at all. 

Still, the success rate for the criterion D is only marginal, mainly because it is intrin- 
sically difficult to jump from -Psat/-fju P = 1-5 to >2.3, and end up below 2.5. First, as 



Morbidelli et al. (2007) was stable. 
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the jump must be rather large, it is hard find cases where Psat/-Pj up is larger than 2.3 and 
smaller than 2.5 after the jump. This is statistically unlikely (3.3% averaging over the first 
seven lines of Table 6) given the spread of possible jump amplitudes for different encounter 
geometries. Second, there is a problem with the residual migration if the planetesimal disk 
is initially massive. In such cases, even if 2.3 < -Psat/-Pju P < 2.5 after the jump, there 
is still a significant mass in the planetesimal disk that needs to be processed through the 
Jupiter /Saturn zone. Consequently, -Psat/-Pj up keeps increasing after the jump and reaches 
> 2.5 (for M disk > 50 M Earth ; e.g., Fig. [7b). The residual migration is difficult to avoid 
unless the planetesimal disk had a relatively low mass to start with. 

An additional problem with the case discussed here is that ess that we obtained in the 
simulations was nearly always too small. This happened because even if ej up was kicked 
by the ejected planet, it was quickly damped at later times by secular friction from the 
planetesimal disk (see, e.g., Fig. [7b). We tested several possible solutions to this problem. 
Disks with M^sk < 35 MEarth did not work for the (3:2,3:2,4:3,5:4) resonant chain because 
two or more planets were lost too often, end even if four planets survived in the end their 
orbits were too wild (B satisfied in 5-10% of cases only). 

We tested dynamically excited planetesimal disks. The disk excitation reduces the disk's 
ability to damp ej up and may thus lead to better results. The disk could have been excited 
by large bodies that formed in it, as evidenced by Pluto-sized bodies in the present Kuiper 
belt. In addition, the disk could have been excited near its inner edge by orbital resonances 
with the outer ice giants. In both cases, the disk planetesimals were distributed according to 
the Rayleigh distribution. We used (e) = 0.15 and (i) = 0.075 (case El), and (e) = 0.1 and 
(i) = (case E2), where brackets denote the mean values. Case El mimicked the excitation 
expected from the large bodies in the disk. Case E2 would be more appropriate for the 
orbital resonances that do not generally excite inclinations. 

Recall that for the cold planetesimal disk we used a minimal separation of A > 0.5 
AU between the outer planet's orbit and semimajor axis of disk particles (see discussion in 
Sect. 2). Now that we deal with excited disks with e up to 0.15, we used A > 3.5 AU for 
consistency, so that the minimal physical distance between planetesimals and planets was 
larger than 0.5 AU. 

With Mdi S k = 50 Msarth, case El gave A = 30% and B = 10%, which was comparable 
to what we obtained for this setup with a dynamically cold disk. Interestingly, there was 
one simulation (out of 30) that matched C as well (e 55 = 0.032). Case E2 gave A = 26% 
and B = 7%, and also one simulation where the criterion C was satisfied (ess — 0.026). 
Conversely, the excited disks with Mdisk = 35 MEarth did not show any case where C would 
be satisfied. We conclude that the disk excitation may help, but the simulations matching 
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C are still disturbingly rare. 

We studied the dependence of the results on the mass of the fifth planet. We found that 
decreasing the fifth planet's mass helps to boost the success rate for A and B. For example, 
with (3:2,3:2,4:3,5:4), Mdisk = 50 MEarth, -8(1), and the mass of the inner ice giant Mi ce = 0.5 
MNep, the success rate in matching A and B was 63% and 30%, while it was 37% and 23% 
with Mi ce = M Nep . In all studied cases with M Ice = 0.5 M Nep , the small inner planet was 
removed (83% ejected and 17% collided with Jupiter). 

On the downside, we found that the jump of Psat/-Pj up and excitation of ej up with 
Mice < 1 MNep were smaller than with Mi ce ~ M^p, which compromised the success rate 
for C and D. In fact, ej up was excited in most of the runs with Mi ce = 0.5 M^ep not by the 
scattering encounter itself, but rather when Jupiter and Saturn crossed the 2:1 resonance. 
As we discussed Sect. 5.1, the resonant crossing has only a minor effect on 655. 

Increasing the mass of the inner ice giant did not help either, because the system became 
violently chaotic during the instability and typically lost two or more planets. In addition, in 
many cases where only one planet was ejected, the lost planet was one of the outer ice giants. 
For example, with = 50 MEarth, #(1) and M\ ce = 2M^ ep we found that A = 15% and 

i? = 3% (when only the removal of the massive inner planet was considered). We conclude 
that the problems with the low success rate of C and D cannot be resolved by changing the 
mass of the fifth planet. 

The problems discussed above may be related to the fact that the ice giants' orbits were 
closely packed together for (3:2,3:2,4:3,5:4), perhaps too closely. We found that the other 
compact resonant chains investigated here, such as (3:2,3:2,4:3,4:3) or even (3:2,3:2,3:2,4:3), 
behave in very much the same way as (3:2,3:2,4:3,5:4), giving way to violent instabilities 
that fail to satisfy the constraints. We therefore decided to focus on more relaxed resonant 
chains. We discuss the results obtained with (3:2,3:2,3:2,3:2) first. These results were more 
encouraging. 

5.2.2. Relaxed Chains with the 3:2 Jupiter- Saturn Resonance 

The (3:2,3:2,3:2,3:2) resonant chain leads to a different mode of instability, if A < 2 AU. 
Figure [S] illustrates a case with M disk = 35 M Earth , 13(1) and A = 1.5 AU. Unlike in all cases 
discussed so far, breaking the resonant lock of the inner ice giant did not result in an im- 
mediate instability. Instead, Neptune migrated deep into the planetesimal disk before any 
scattering encounters between planets occurred. As Neptune and Uranus scattered planetes- 
imals into the Jupiter/Saturn zone, Jupiter and Saturn underwent divergent migration and 
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left the 3:2 resonance. Their eccentricities were then excited by the the 5:3 resonance cross- 
ing, and the instability happened shortly after (16.8 Myr after the start of the simulation). 

The epoch of planetary encounters was brief, lasted from 16.816 Myr to 16.9 Myr 
at which point the fifth planet was ejected from the Solar System (escape speed = 
0.74 km s-MQ All planets participated in the encounters, as required for capture of the 
irregular satellitesjfl but all these encounters involved the fifth planet only. None of the sur- 
viving giant planets had an encounter with another surviving giant. In total, 282 encounters 
happened in this case, where an encounter is defined here by the condition that the Hill 
spheres of the two planets at least partially overlap. 

Now, the case illustrated in Fig. simultaneously matched all our criteria. Amplitude 
e 55 was excited by the ejection of the fifth planet by Jupiter and was not damped too much 
during the following evolution, because the planetesimal disk had a relatively low initial mass, 
and was partially disrupted by Neptune prior to the excitation event. The final amplitudes 
were = 0.024 and = 0.012. Here, the amplitude was lower than real one, but in 
another simulation (with the same setup) that matched all our criteria as well, we obtained 
e 55 = 0.041 and e 56 = 0.013. 

The scattering encounters with the fifth planet produced a jump of Psat/-Pjup from 1.7 
to 2.5, just as required to avoid the secular resonances with the terrestrial planets (Fig. 
Efc). There was a small problem with the residual migration in this run, because -Psat/-Pjup 
continued to evolve past 2.5, while Psat/^Jup = 2.49 in the present Solar System, but we 
obtained better results in another simulation with the same setup (Figure |9j although in this 
case was damped to 0.013). 

The overall distribution of planetary orbits obtained in simulations with the (3:2,3:2,3:2,3:2) 
resonant chain, Mdi S k = 35 Mgarth, S(l) and A = 1.5 AU is shown in Fig. [TUJ The overall 
success rate in matching the criteria A, B, C and D was 33%, 16%, 4% and 8%, respectively 
(note that matching C or D required that B was matched). The overall success rate for 
simultaneous matching of C and D was 4%. This is the best result discussed so far. For a 
comparison, a similar setup with a more massive disk, M disk = 50 M Earth , gave to A = 30%, 
B = 17%, C = 0% and D = 3%. The success for C dropped because e 55 was damped more 
efficiently with the massive planetesimal disk. A light disk with Mdi S k = 20 Msarth gave 
A = 20%, B = 7%, C = 3% and D = 0%. 



7 Typically, = 0.5-2 km s" 1 . 

8 We do not define planetary encounters as a separate criterion, because constraints C and D require 
planetary encounters and are generally much more restrictive. 



-20 - 



The mode of instability discussed for (3:2,3:2,3:2,3:2) above may be problematic because 
the migration of Uranus and Neptune prior to the instability may require that the inner edge 
of the planetesimal disk was initially close to the outer ice giant. The edge may then need 
to be fine-tuned to generate the delay between the formation of the Solar System and the 
LHB. This is a problem similar to that of the original Nice model, which prompted Levison 
et al. (2011) to consider distant disks and a different trigger mechanism. Here we do not 
study the problem of delay in detail. Instead, given that the instability mode with prior 
migration of Neptune works much better than any other case, we pursued investigations into 
this instability mode further. 

We first tested several cases where the inner ice giant was placed the 2:1 resonance 
with Saturn ((3:2,2:1,2:1,3:2), (3:2,2:1,3:2,2:1), etc.). We found that these resonant chains 
did not work, because the orbits spread without suffering any major instability, which left 
five planets behind. Unlocking the resonance of one of the ice giants did not destabilize the 
system as well, mainly because the inner ice giant was initially far from Saturn, and safe, 
even if its orbit was not locked deep in the 2:1 resonance. 

Motivated by these results, we considered several resonant chains with the inner ice 
giant in the exterior 3:2 (or 4:3) resonance with Saturn, and in the 2:1 resonance with the 
inner surviving ice giant (e.g., (3:2,3:2,2:1,2:1), (3:2,4:3,2:1,2:1) and (3:2,3:2,2:1,3:2). This 
initial setup increases the likelihood that the inner ice giant gets ejected, because it starts 
relatively close to Saturn, and that the outer two ice giants survive, because there is initially 
a large radial gap between them and the inner ice giant. 

We found that the extended resonant chains such as (3:2,3:2,2:1,2:1) and (3:2,4:3,2:1,2:1) 
did not work because Neptune migrated beyond 30 AU, independently of the assumed mass 
and extension of the planetesimal disk. The explanation for this is clear. As in all simu- 
lations, Neptune scattered planetesimals outward and inward relative to its orbit. Many of 
those that were scattered outward were scattered to a > 30 AU, because initial a^ep — 25 
AU for the extended resonant chains. The ones that were scattered inward must have been 
scattered relatively strongly to be handed to Uranus, because the orbits of Uranus and Nep- 
tune were radially spaced (due to the initial 2:1 resonance). As a result, Neptune migrated 
outward, being propelled by the strong scattering encounters, all the way into the cloud of 
planetesimals that it scattered beyond 30 AU. This occurred even if the planetesimal disk 
had low mass (Mdi S k = 20 Msarth) and/or was initially narrow (spanning, e.g., 26-27 AU). 
Also, even lighter disks with M disk < 15 M Earth did not show much success for A. 

Figure [TT] shows the distribution of final orbits for the (3:2,3:2,2:1,3:2) resonant chain 
and Mdisk = 20 Msarth- These results are more promising. As the distributions obtained with 
B(0) and 13(1) did not differ (assuming A ~ 1 AU), we combined them together in Fig. [TTJ 
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The combined success rate was A = 32%, B = 12%, C = 7%, D = 5% and C*&£> = 5%. This 
was similar to the statistics obtained for (3:2,3:2,3:2,3:2) and M disk = 35 M Earth (Fig. [TO]) . 
but here things were slightly better in detail. 

For example, the distribution of model aura nicely matched the present semimajor axis 
of Uranus, while the model aura was slightly lower than required in Fig. [TUJ The model values 
of ej up obtained with (3:2,3:2,2:1,3:2) were almost perfect in that they closely surrounded 
the target ej up = 0.046. Figures [T2"| [TBI and IT41 show examples of the successful simulations. 
We illustrate these cases in detail because they were some of the best results obtained in 
this work. 

The simulation in Fig. [U ((3:2,3:2,2:1,3:2) and M disk = 20 M Earth ) ended with g 5 = 4.95 
arcsec yr" 1 , g 6 = 29.7 arcsec yr _1 , g 7 = 4.53 arcsec yr _1 , and g s = 0.52 arcsec yr _1 . These 
values were a close match to those listed for the real planets in Table 2. The amplitudes 
of eccentricity modes were = 0.049, = 0.033, = 0.024, = 0.041, = 0.042, 
e67 = 0.003, and = 0.039. The frequencies and amplitudes of inclination modes were 
equally good. Neptune ended with an eccentricity that was about twice as large as its 
present value, but this was at least partly related to a minor mismatch in the simulation, 
because Uranus and Neptune crossed the 2:1 resonance, which did not happened in reality. 

Figure [13] shows another case obtained with the same setup. In this case, g$ = 4.83 
arcsec yr _1 , g§ = 28.9 arcsec yr -1 , g 7 = 3.55 arcsec yr _1 , and gs = 0.64 arcsec yr _1 . The 
amplitudes were e 55 = 0.027, e 56 = 0.014, e 57 = 0.002, e 65 = 0.041, e 66 = 0.024, e 67 = 0.002, 
677 = 0.033. Finally, the case illustrated in Fig. [lj] had frequencies g$ = 4.66 arcsec 
yr _1 , g Q = 28.1 arcsec yr _1 , g-j = 3.45 arcsec yr _1 , g$ = 0.61 arcsec yr" 1 , and amplitudes 
e 55 = 0.024, e 56 = 0.019, e 57 = 0.004, e 65 = 0.019, e 66 = 0.057, e 67 = 0.004, e 77 = 0.040, 
eg8 = 0.007. This case differs with respect to the previous two in that the innermost ice 
giant survived and evolved onto Uranus-like orbit, while the middle ice giant was ejected. 
The Neptune's model eccentricity was very good in this simulation. 

Using (3:2,3:2,2:1,3:2) and Mdi S k — 20 MEarth we tested how the results were affected by 
small changes (up to 50%) in the mass of the fifth planet. We found that it was probably not 
lower than 0.7 M Ura because in the simulations with these low masses Jupiter (and Saturn) 
were not kicked enough by the planet's ejection. Masses slightly larger than MN ep showed 
lower success rates for A and B, but still worked relatively well for C and D. To summarize, 
we find that the fifth planet probably had mass in the Uranus/Neptune range, which is 
encouraging in that we do not need to invoke any special mass regime. 
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5.2.3. 2:1 Jupiter- Saturn Resonance 

The initial conditions discussed in this section are probably academic, because the 2:1 
resonance between Jupiter and Saturn is not favored by the standard model of the migration 
of planets in the protoplanetary gas disks, and their capture in resonances (see Sect. 3). We 
therefore discuss this case briefly, concentrating on the main differences with respect to the 
results obtained with more conservative assumptions. 

The five-planet case with the 2:1 resonance would require that the (ejected) inner ice 
giant had lower mass, because when Jupiter and Saturn start in the 2:1 resonance, their 
period ratio needs to change by ~0.5 only (from 2 to 2.49), which requires a smaller pertur- 
bation. The best results were obtained with M\ ce = 0.5 MN ep and 15 < Mdi s k < 35 MEarth- 
For example, with Mdisk = 20 MEarth, we obtained a ~50-70% success for A and 20-40% 
success for B. In addition, because the required jump of -Psat/-fju P is smaller and easier to 
achieve, the simulations also show a large success for D (reaching 20% in the most favorable 
cases) . 

Figure [TBI illustrates this case. The model semimajor axis of Uranus and Neptune were 
a bit smaller than what would be ideal, but this should easily be corrected by extending the 
disk slightly beyond 30 AU (we used r out = 30 AU). All else looked great, except that the 
Jupiter's eccentricity was generally too small and violated constraint C. We found only two 
simulations out of more than 1,000, where the criterion C was satisfied. 

This problem arises because while a relatively low-mass ice giant is required to produce 
the needed small jump of -Psat/-Pjup ; the same low-mass planet is generally incapable of ex- 
citing e 55 enough during encounters. In addition, the problem with secular friction discussed 
in the previous sections is also apparent here. To show this clearly, we adjusted the initial 
system so that Jupiter had a substantial eccentricity initially (Fig. [16]). This did not helped 
at all because the initial eccentricity was quickly damped. 

The problem with matching C in the case with five planets and 2:1 resonance between 
Jupiter and Saturn is difficult to avoid (at least we were not able to resolve this problem 
with 1,000+ simulations). Thus, even if the success rate for A, B and D was promisingly 
large, we are not overly optimistic about this case. A more thorough testing will be required 
to reach a more definitive conclusion!^ 



9 The discussion presented here was based on the overall synthesis of our simulation results. Note that it 
can be misleading to compare two specific simulation sets in Table 6. For example, we have just been lucky 
to find one case (out of 30) for (2:1,3:2,3:2,3:2), Mdisk = 20 MEarth and B(l) for which the criterion C was 
satisfied. In contrast, the results for the 3:2 Jupiter- Saturn resonance were consistently good, showing >1% 
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5.3. Case with 6 Planets 

We performed 1290 simulations in total for six initial planets and nine different resonant 
chains. Given that the number of parameters in the six-planet case is much larger than 
the one with four or five planets, we are not confident that we sampled parameter space 
exhaustively enough to make detailed conclusions. Still, the six-planet case is very interesting 
because we were able to obtain good results even after a relatively small number of trials. 

Figure [17] illustrates one of the most promising six-planet results. This result was 
obtained for the (3:2,4:3,3:2,3:2,3:2) resonant chain and Mdi S k = 20 Mgarth- The two extra 
ice giants were placed between Saturn and the inner surviving ice giant, and were given 
masses M\ ce i = Mi ce2 = 0.5 MN ep (the six-planet case with Mi ce i = Mi ce2 = 1 MN ep does not 
work, because the instability is too violent). 

The results show a large spread, because the six planets have complex interactions 
during the instability. The distributions of model orbits nicely overlap with the real orbits. 
Even the details match. For example, most simulations ended with eN — 0.01, just as needed. 
Given our previous experience with the simulations of instability in the four- and five-planet 
cases, the agreement in Fig. [T7J is remarkable. Because of the larger spread, however, the 
success rate was lower than that for our best five-planet simulations. In the specific case 
discussed above, we obtained A = 30%, B = 10%, C = 3%, D = 3% and CkD = 2%. 

The instability typically occurred in two steps, corresponding to the ejection of the two 
planets. Sometimes, as in Fig. [T8J the ejection of the two planets was nearly simultaneous, 
but most of the times there was a significant delay between ejections. This was useful because 
the first planet's ejection partially disrupted the planetesimal disk and reduced its capability 
to damp ess, which was then excited by the second planet's ejection. While this mode of 
instability can be important, we would need to increase the statistics (>100 simulations 
for each initial condition) to be able to properly resolve the small success fractions in the 
six-planet case. 



6. Discussion 

One of the main results discussed in this paper is that the Solar System instability with 
five initial planets tends to give better results then the instability with four initial planets. 
Batygin et al. (2012; hereafter B12) performed a similar analysis and found instead that the 



success rate for matching all criteria simultaneously. 
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four- and five-planet cases show about the same success rate in matching constraints. Here 
we discuss the possible causes of this disagreement. 

B12 used an iV-body integrator with forces that mimic the effects of gas to generate 
the initial resonant chains of planets. This method is similar to that described for our 
Phase-1 integrations in Sect. 2, and should result in similar initial orbits for Phase 2. The 
resonant chains explicitly discussed in B12 included (3:2,2:1,5:4,5:4) and (3:2,3:2,2:1,4:3). 
We discussed this type of relaxed resonant chains in Sect. 5.2.2. The initial properties of 
the planetesimal disk in B12 were also similar to the ones used here. We therefore believe 
that the initial conditions are not the main cause of the disagreement. Instead, we explain 
below that different conclusions were reached in B12 probably because B12 gave a different 
emphasis to different constraints (see Sect. 3 for our constraints). 

On one hand, B12 found, correcting the results of Batygin & Brown (2010), that ~ 10% 
of simulations with four initial planets matched constraints. Their constraints, however, 
were somewhat different and less restrictive than the ones we used here. First of all, B12 
did not apply any upper limit on the excitation of planetary orbits, while we required in the 
criterion B that e < 0.11 and % < 2° (relative to the invariant plane). It is therefore possible 
that in at least some of the simulations of B12, which gave reasonable e^, some planetary 
orbits had e > 0.11 and/or i > 2°. 

B12 did not consider the constraint from the terrestrial planets (our criterion D). It is 
therefore possible that in some of the simulations that were found to be successful in B12, 
the evolution of the mode was smooth, which would lead to the secular resonances with 
the terrestrial planets. Note that some of these models could still be valid, but the fraction 
should be small, probably <10% (Brasser et al. 2009). 

While the two issues pointed out above can resolve some of the main differences between 
B12 and our results, we are still puzzled, given our clearly negative results for the four-planet 
case, that B12 were able to find positive results for the four-planet case even with very massive 
disks (Mdi S k > 50 Msarth)- We tried to repeat the simulations reported in B12 and found that 
massive disks efficiently damp and lead to problems with residual migration of Jupiter 
and Saturn. 

On the other hand, B12 considered constraints from the Kuiper belt that we ignored 
here. As B12 pointed out, only about half of their five-planet simulations were compatible 
with the low eccentricities and low inclinations in the cold classical Kuiper belto while 



The cold classical Kuiper Belt is a population of trans-Neptunian bodies dynamically defined as having 
orbits with semimajor axis a = 42-48 AU, perihelion distances that are large enough to avoid close encounters 
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most of their four-planet cases were fine. We may have been therefore overly optimistic, by 
a factor of ~2, in favoring the five-planet case over the four-planet case. The constraints 
on the planetary instability from the small body reservoirs in the Solar system (Kuiper belt 
objects, asteroids, Trojans and satellites) are clearly very important. We will consider these 
constraints in the follow-up work. 

B12 and this work agree on one central issue. As B12 was not concerned with the 
delay between the protoplanetary nebula dispersal and LHB, they placed the planetesimal 
disk's inner edge just beyond the outer ice giant's orbit. This triggered, almost immediately, 
Neptune's fast migration through the disk and the instability mode akin to that we discussed 
for our most successful runs in Sect. 5. B12 therefore incidentally explored the setup that 
we favored here based on a broader sampling of the initial parameters. It remains to be 
understood whether this setup can lead to the late instability without overly restrictive 
assumptions on the structure of the planetesimal disk. 

It would be useful in this context, for example, if the LHB started ~ 4.2 Gyr ago as sug- 
gested by Bottke et al. (2012), because a shorter time delay since the protoplanetary nebula 
dispersal (~350 Myr) would help to relax constraints on A. Alternatively, the planetesimal 
disk could have been stirred by large bodies that formed in it, and spread over time, so that 
the effective A decreased. This could lead to the late instability even if the initial A was 
large. Investigations of these issues are left for future work. 

7. Conclusions 

Recent studies suggest that Jupiter and Saturn formed and migrated in the protoplane- 
tary gas disk to reach a mutual resonance, most likely the 3:2 resonance, where the Saturn's 
orbital period was 3/2 longer than that of Jupiter. After the gas disk's dispersal, the orbits 
of Jupiter and Saturn must have evolved in some way to eventually arrive to the current 
orbits with the orbital period ratio of 2.49. This can most easily be achieved, considering 
constraints from the terrestrial planets and e^, if Jupiter and Saturn scattered off of Uranus 
or Neptune, or a planet with mass similar to that of Uranus or Neptune. 

We performed iV-body integrations of the scattering phase between the Solar System's 
giant planets, including cases where one or two extra ice giants were assumed to have formed 
in the outer Solar System and ejected into interstellar space during instability. We found 
that the initially compact resonant configurations and low masses of the planetesimal disk 



to Neptune, and low inclinations (i < 5°). 
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(M disk < 50 M Ea rth) typically lead to violent instabilities and planet ejection. On the other 
hand, the initial states with orbits that are more radially spread (e.g., Jupiter and Saturn in 
the 2:1 resonance) and larger Mdi s k result in smooth migration of the planetary orbits that 
leads to incorrectly low and excitation of the terrestrial planet orbits. Finding the sweet 
spot between these two extremes is difficult. 

Some of the statistically best results were obtained when assuming that the Solar System 
initially had five giant planets and one ice giant, with the mass comparable to that of Uranus 
and Neptune, was ejected to interstellar space by Jupiter. The best results were obtained 
when the fifth planet was assumed to have the mass similar to Uranus/Neptune, was placed 
on an orbit just exterior to Saturn's (3:2 and 4:3 resonances work best), and the orbits 
of Uranus and Neptune migrated into the planetesimal disk before the onset of planetary 
scattering. This mode of instability is favored for several reasons, as described below. 

As planetesimals are scattered by Uranus and Neptune and evolve into the Jupiter/Saturn 
region, Jupiter, Saturn and the fifth planet undergo divergent migration. This triggers an 
instability during which the fifth planet suffers close encounters with all planets and is even- 
tually ejected from the Solar System by Jupiter. Uranus and Neptune generally survive 
the scattering phase, because their orbits migrated outward during the previous stage and 
opened a protective gap between them and the gas giants. This mode of instability produces 
just the right kind of Jupiter's semimajor axis evolution -known as jumping Jupiter- that is 
required from the terrestrial planet constraint. 

Moreover, e^, excited by the fifth planet ejection, is not damped to incorrectly low 
values by secular friction from the planetesimal disk, because the planetesimal disk had 
been disrupted by Uranus and Neptune before the excitation event. The low mass of the 
planetesimal disk at the time of planet scattering also leads to only a brief migration phase of 
Jupiter and Saturn after the scattering phase, and prevents -Psat/-Pjup from evolving beyond 
its current value. The excessive residual migration of Jupiter and Saturn was a problem in 
most other cases investigated here. 

The mode of instability with early migration of Uranus and Neptune is problematic, 
however, because the inner edge of the planetesimal disk may need to be fine-tuned to 
generate the delay between the formation of the Solar System and the LHB. In addition, the 
range of possible outcomes is rather broad, indicating that the present Solar System is neither 
typical nor an expected result, and occurs, in best cases, with only a ~5% probability (as 
defined by simultaneous matching of all four criteria defined in Sect. 3). In ~95% of cases, 
the simulations ended up failing at least one of our constraints. This may seem unsatisfactory, 
but given the issues discussed in Sect. 4, the fact that our criteria are relatively strict, and 
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because the instability-free model does not work at alio these findings should be seen in a 
positive light. An important follow-up of this work will be to consider additional constraints 
from the small body populations (e.g., the dynamical structure of the Kuiper belt), and see 
whether the successful simulations identified here will also match those constraints. 

The case with six giant planets is also interesting in that the instability occurs in two 
steps, corresponding to the ejection of the two planets. The best results were obtained in 
this case when the two ejected planets were given similar masses (about half the mass of 
Neptune) and were placed between the orbits of Saturn and the inner surviving ice giant. 
As expected, the six-planet case leads to a larger variety of results than the five-planet 
case. The probability of ending the six-planet simulation with the present properties of the 
solar system is therefore lower than in the five-planet case. Still, our six-planet results are 
fundamentally better than those obtained in the four-planet case, where the differences were 
systematic (e.g., never large enough), and the success rate was below the resolution limit 
of our study. 
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a (AU) 


e 




Jupiter 


5.203 


0.046 


0.37 


Saturn 


9.555 


0.054 


0.90 


Uranus 


19.22 


0.044 


1.02 


Neptune 


30.11 


0.010 


0.67 



Table 1: Mean orbital elements of planets. The mean values reported here were obtained 
by numerically integrating the orbits of Jupiter, Saturn, Uranus and Neptune for 10 Myr, 
and computing the average of orbital elements over this interval. The mean inclinations are 
given with respect to the invariant plane. 



j Qj (arcsec yr v ) Sj (arcsec yr x ) 

5 4.24 

6 28.22 -26.34 

7 3.08 -2.99 

8 0.67 -0.69 



Table 2: Secular frequencies of giant planets in the Solar System. The frequencies were 
obtained by numerically integrating the orbits for 10 Myr, and Fourier analyzing the results. 
The S5 frequency vanishes when the inclinations are referred to the invariant plane. 



5 6 7 8 

Jupiter 0.044 0.015 0.002 

Saturn 0.033 0.048 0.002 

Uranus 0.038 0.002 0.029 0.002 

Neptune 0.002 - 0.004 0.009 



Table 3: Secular amplitudes e^, where j denotes individual planets (5 to 8 from Jupiter to 
Neptune). The proper modes of each planet's orbit are denoted in bold. The amplitudes 
were obtained by numerically integrating the orbits for 10 Myr, and Fourier analyzing the 
results. Values lower than 0.001 are not reported. 
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6 


7 


8 


Jupiter 


0.36 


0.06 


0.07 


Saturn 


0.90 


0.05 


0.06 


Uranus 


0.04 


1.02 


0.06 


Neptune 




0.12 


0.66 



Table 4: Secular amplitudes ijk- The proper modes of each planet's orbit are denoted in bold. 
The amplitudes are reported in degrees. They were obtained by numerically integrating the 
orbits for 10 Myr, and Fourier analyzing the results. Values lower than 0.01° are not reported. 
The S5 frequency vanishes in the invariant plane and gives no contribution here. 
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Table 5: The results of selected four-planet models. The columns are: the (1) mass of the 
planetesimal disk (Mdisk), (2) A defined as r in — a n , where r in is the initial radial distance of 
the inner edge of the planetesimal disk and a n is the semimajor axis of the outer ice giant, 
(3) radial distance of the outer edge of the planetesimal disk (r out ), (4) B(j) specifying the 
instability trigger (see Sect. 2 for a definition), (5) number of simulations performed for each 
case (Ajob), and (6-9) success rate for our four criteria defined in Sect. 3. 



- 34 - 



Miisk 


A r out B(j) 


Nsim A 


B 


C 


D 


(M E arth) 


(AU) (AU) 


% 


% 


% 


% 




(3:2,3:2,4:3,5:4), a 5 


= 13.9 AU 








35 


1.0 30 1 


30 13 


3 





3 


50 


1.0 30 1 


30 37 


23 








50 


3.5 30 1 


30 23 


10 





3 




(3:2,3:2,3:2,3:2), a 5 


= 16.1 AU 








35 


1.5 30 


30 23 


3 





3 


35 


1.5 30 1 


100 33 


16 


4 


8 


50 


1.5 30 1 


100 30 


17 





3 




(3:2,3:2,4:3,4:3), a 5 


= 14.5 AU 








50 


1.0 30 1 


30 47 


23 


3 


3 




(3:2,3:2,2:1,3:2), a 5 


= 22.2 AU 








20 


1.0 30 


30 33 


13 


7 


3 


20 


1.0 30 1 


30 30 


10 


7 


7 


35 


1.0 30 1 


30 33 


17 





10 




/000001 O 1 \ 

(3:2,3:2,2:1,2:1), a 5 


O A r ATT 

= 24.5 AU 








35 


1.0 30 


30 43 


17 


7 


7 


35 


1.0 30 1 


30 23 


13 


3 


3 


35 


2.0 30 1 


30 30 


23 


3 


3 


35 


3.0 30 1 


100 44 


19 


3 


3 




(2:1,3:2,3:2,3:2), a 5 


= 19.3 AU 








20 


1.0 30 1 


30 53 


36 


3 


20 


35 


1.0 30 1 


30 53 


43 





17 




(2:1,4:3,3:2,3:2), a 5 


= 17.9 AU 








20 


1.0 30 1 


30 67 


42 





17 


20 


3.5 30 1 


30 75 


25 





20 



Table 6: The results of selected five-planet models. See the caption of Table 5 for a definition 
of different parameters shown here. 
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Table 7: The results of selected six-planet models. See the caption of Table 5 for a definition 
of different parameters shown here. 
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Fig. 1. — Orbit histories of the giant planets in a simulation with four initial planets. The 
four planets were started in the (3:2,4:3,3:2) resonant chain, M disk = 35 M Earth and B(l) (see 
Sect. 2 for the definition of B(l)). (a) The semimajor axes (solid lines), and perihelion and 
aphelion distances (dashed lines) of each planet's orbit. The red, green, turquoise and blue 
lines correspond to Jupiter, Saturn, Uranus and Neptune. The black dashed lines show the 
semimajor axes of planets in the present Solar System, (c) The period ratio -Psat/-Pj up - The 
dashed line shows -Ps a t/-Pju P = 2.49, corresponding to the period ratio in the present Solar 
System. The shaded area approximately denotes the zone where the secular resonances with 
the terrestrial planets occur, (b) Jupiter's eccentricity, (d) Jupiter's inclination. The dashed 
lines in (b) and (d) show the present mean eccentricity and inclination of Jupiter. 
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Fig. 2. — Final orbits obtained in our simulations with four planets started in the (3:2,3:2,4:3) 
resonant chain, M disk = 75 M Earth and B(l) (see Sect. 2 for the definition of £>(1)). (a) Mean 
eccentricity, (b) Mean inclination. The mean orbital elements were obtained by averaging 
the osculating orbital elements over the last 10 Myr (i.e., from 90 to 100 Myr). Only the 
systems ending with four planets are plotted here (dots). The bars show the mean and 
standard deviation of the model distribution of orbital elements. The mean orbits of real 
planets are shown by triangles. Colors red, green, turquoise and blue correspond to Jupiter, 
Saturn, Uranus and Neptune. 
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Fig. 3. — Orbit histories of the giant planets in a simulation with four initial planets. See 
the caption of Fig. [TJfor the description of orbital parameters shown here. The four planets 
were started in the (3:2,3:2,4:3) resonant chain, and Mdisk = 75 MEarth- 
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Fig. 4. — Orbit histories of the giant planets in a simulation with four initial planets. See 
the caption of Fig. [TJfor the description of orbital parameters shown here. The four planets 
were started in the (3:2,3:2,4:3) resonant chain, M^isk = 50 MEarth and A = 1 AU. 
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Fig. 5. — Final orbits obtained in our simulations with four planets started in the (2:1,3:2,3:2) 
resonant chain, and Mdisk = 35 Msarth- See the caption of Fig. [2] for the description of orbital 
parameters shown here. 
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Fig. 6. — Final orbits obtained in our simulations with five planets started in the 
(3:2,3:2,4:3,5:4) resonant chain, Mdisk = 50 Msarth and A = 1 AU. See the caption of Fig. [2] 
for the description of orbital parameters shown here. 
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Fig. 7. — Orbit histories of the giant planets in a simulation with five initial planets. See 
the caption of Fig. Q] for the description of orbital parameters shown here. The five planets 
were started in the (3:2,3:2,4:3,5:4) resonant chain, M disk = 50 M Earth and A = 1 AU. The 
fifth planet was ejected at t = 0.8 Myr after the start of the simulation. 
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Fig. 8. — Orbit histories of the giant planets in a simulation with five initial planets. See 
the caption of Fig. Q] for the description of orbital parameters shown here. The five planets 
were started in the (3:2,3:2,3:2,3:2) resonant chain, M disk = 35 M Earth and 13(1). The fifth 
planet was ejected at t = 17 Myr after the start of the simulation. 



-44- 




Time (yr) 



to 



05 
o 



S CM 



u 

4> 

a, 



1 


1 — 1 — 1 1 1 1 1 1 1 1 — 1 — 1 1 1 1 1 1 


— 1 1 — 1 1 1 1 1 1 

(c) : 


1 1 1 












1 


1 1 1 1 1 1 * 





10° 



10° 



10' 



10° 



a 

« d 

4) ° 

O 

O 

w 



1 1 — 1 — 1 1 1 1 1 1 1 1 — 1 — 1 1 1 1 1 1 


- 1 1 — 1 — 1 1 1 1 1 

(b) ■ 


• 





10° 



10° 



10' 



10° 



Time (yr) 



^ 00 

so d 
« 

' CO 

13 d 
o 

■l-j 

* * 

a o 
■i-i 

"o 

- d 



1 1 — 1 — 1 1 1 1 1 1 1 1 — 1 — 1 1 1 1 1 1 


— 1 1 — 1 1 1 1 1 1 

(d) | 








10= 



10° 



10' 



10° 



Time (yr) 



Time (yr) 



Fig. 9. — Orbit histories of the giant planets in a simulation with five initial planets. See the 
caption of Fig. [T]for the description of orbital parameters shown here. The five planets were 
started in the (3:2,3:2,3:2,3:2) resonant chain, M disk = 35 M Earth and 25(0) The fifth planet 
was ejected at t = 15.2 Myr after the start of the simulation. 



-45 - 




-i — i — i — i — | — i — i — i — i — | — i — i — i — i — | — i — i- 

(b) 



tuO 
03 

"ti CM - 

Hi 
O 

• rH 

cd 

• rH . — 
i 1 

o 

rH 



A 



O 







_i i i i I i i i i I i i i l_ 



A . 



J I u 



10 20 
Semimajor axis (AU) 



30 



Fig. 10. — Final orbits obtained in our simulations with five planets started in the 
(3:2,3:2,3:2,3:2) resonant chain, M disk = 35 M Earth and 13(1). See the caption of Fig. [2] 
for the description of orbital parameters shown here. 



- 46 - 



-i 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 r 



O 

h-> 
• i— i 

O 

•i-H 
S-H 

-M 

CD O - 

o 
o 
H 



o 



IP 



(9) 



o 







10 



20 



in 
o 

•i-H 

■+J 

cd 

• i-H , — 

a 



A 

I 

% 



A 



_l I I I I I I I I I I I I I I 1_ 



30 



~ i 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 r 



(b) 



A 



o 



j i i i i i i i i i i i i i i i i_ 







10 20 
Semimajor axis (AU) 



30 



Fig. 11. — Final orbits obtained in our simulations with five planets started in the 
(3:2,3:2,2:1,3:2) resonant chain, and M disk = 20 M Earth . Results obtained with B(0) and 
B(l) for A = 1 AU were combined here. See the caption of Fig. [2] for the description of 
orbital parameters shown here. 
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Fig. 12. — Orbit histories of the giant planets in a simulation with five initial planets. See 
the caption of Fig. Q] for the description of orbital parameters shown here. The five planets 
were started in the (3:2,3:2,2:1,3:2) resonant chain, M disk = 20 M Earth and 13(1). The fifth 
planet was ejected at t = 34.4 Myr after the start of simulation. 
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Fig. 13. — Orbit histories of the giant planets in a simulation with five initial planets. See 
the caption of Fig. Q] for the description of orbital parameters shown here. The five planets 
were started in the (3:2,3:2,2:1,3:2) resonant chain, M disk = 20 M Earth and B(0). The fifth 
planet was ejected at t = 14.9 Myr after the start of the simulation. 
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Fig. 14. — Orbit histories of the giant planets in a simulation with five initial planets. See 
the caption of Fig. Q] for the description of orbital parameters shown here. The five planets 
were started in the (3:2,3:2,2:1,3:2) resonant chain, M disk = 20 M Earth and 13(1). The fifth 
planet was ejected at t = 6.1 Myr after the start of the simulation. 
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Fig. 15. — Final orbits obtained in our simulations with five planets started in the 
(2:1,3:2,3:2,3:2) resonant chain, and Mdi S k = 20 MEarth- See the caption of Fig. |2] for the 
description of orbital parameters shown here. 
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Fig. 16. — Orbit histories of the giant planets in a simulation with five initial planets. See 
the caption of Fig. Q] for the description of orbital parameters shown here. The five planets 
were started in the (2:1,3:2,3:2,3:2) resonant chain, and Mdi S k = 20 Mgarth- The fifth planet 
was ejected at t — 5.6 Myr after the start of the simulation. 
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Fig. 17. — Final orbits obtained in our simulations with six planets started in the 
(3:2,4:3,3:2,3:2,3:2) resonant chain, and Mdisk = 20 Msarth- See the caption of Fig. |2]for the 
description of orbital parameters shown here. 



-53 - 



o 
-a- 



3 



« 

o 
'5* 

s 

•l-H 

a 

a) 

CO 



CM 



1 1 — i — i i i 1 1 1 r 

1 


n — i i i 1 1 1 1 1 — i — i i i 1 1 










r- 




f - - - ' z " = - - ~* - - tyVXZXSSSsaaammm 
■ inn* : 


i i i * i i i * i i i 



10° 



10° 



10' 



10° 



Time (yr) 



m 

<N 



<d 
OS 

■o 
o 

•c 

a. 



1 


i — i — i i i 1 1 1 1 r 


— i — i i i i 1 1 1 1 — i — i i i i i 

( C ) : 














■ 


1 ■ ■ 


i 



10 3 



10° 



10' 



10° 



T 1 1 I I I I I I 1 1 1 I I I I I I 1 1 1 I I I I I 




- 1 1 — i — i i i 1 1 1 r 



10° 10' 10° 

Time (yr) 

17TTTT1 1 — 1 — 1 — 1 1 1 1 




Time (yr) 



Time (yr) 



Fig. 18. — Orbit histories of the giant planets in a simulation with six initial planets. See 
the caption of Fig. [1] for the description of orbital parameters shown here. The six planets 
were started in the (3:2,4:3,3:2,3:2,3:2) resonant chain, and Mdi S k = 20 Msarth- The fifth and 
sixth planet were ejected at t — 3.18 and 3.65 Myr after the start of simulation (yellow and 
purple lines in (a)). 



